Stochastic resonance in a non Markovian discrete state model for excitable systems 
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We study a non Markovian three state model, subjected to an external periodic signal. This 
model is intended to describe an excitable systems with periodical driving. In the limit of a small 
amplitude of the external signal we derive expressions for the spectral power amplification and the 
signal to noise ratio as well as for the inter-spike interval distribution. 
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Stochastic resonance (SR) is one of the interesting ef- 
fects where noise plays a constructive role in nature. It 
describes the optimal response of a nonlinear system to 
a periodic signal at a finite non zero noise level. While 
today most research in the field of SR is focused on ap- 
plications in information processing systems like neurons, 
the effect was originally discovered in the context of cli- 
matology 0, where it allowed to explain the periodic 
occurrence of ice ages. 

A successful theory of SR in bistable systems was intro- 
duced in • There a reduction to a two state Markovian 
dynamic has been proven to be very useful to obtain an- 
alytical expressions of the signal to noise ratio (SNR), 
which as a hallmark of SR shows a maximum at a cer- 
tain finite noise level. Another example of a two state 
theory for a bistable system with delay has been devel- 
oped in Q. In this system the spectral power amplifi- 
cation (SPA) shows several maxima as a function of the 
frequency. Array enhanced SR in coupled bistable sys- 
tems is studied in Qj]. Recently, in || a theory of SR 
in non Markovian two state models has been established 
and applied to ion channel gating dynamics. In 0] a 
piecewise linear FitzHugh-Nagumo system was modeled 
as a two state dynamics which allowed to calculate the 
SPA in this particular case. 

In this paper study SR in excitable systems [?1 Isl IflL EH 
Hl] | by mapping the characteristic parts of an excitable 
stochastic dynamics onto a non Markovian three state 
dynamics [JjJ]. The resulting model allows to derive for- 
mulas for relevant quantities like SNR or SPA, which may 
help to explain the SR found in real world excitable sys- 
tems driven by periodic signals, ranging from crayfish 
mechanorece pto rs Q to semiconductor lasers with opti- 
cal feedback 

A key feature of excitable stochastic systems is a stable 
fixed point or rest state from which the system can be 
excited by noise onto an excitation loop. This fixed point 
corresponds to state 1 in our model. The transition out 
of this state is assumed to be a rate process, describing 
the excitation over a threshold due to noise. The other 
two states model two parts of the excitation loop each 
having a different output. One might for example think 
of these two parts as the firing and refractory state in 



the dynamics of an excitable FitzHugh-Nagumo neuron 
where the output is high in the firing and low in the re- 
fractory state. The transitions from state 2 to 3 and 3 
back to 1 are governed by arbitrary waiting time distri- 
butions (WTD), which must be adapted to the system 
to be described by this model. 

The external periodical driving acts onto the system 
as a periodical modulation of the rate for the transition 
1 — > 2. The times spent in the firing and refractory 
state are assumed not to be affected by the external driv- 
ing, which is true for example in a stochastic FitzHugh- 
Nagumo system with small driving amplitudes and suffi- 
ciently low driving frequencies. 




Figure 1: Three state model of an excitable system. The 
transition from state 1 to 2 is governed by a Poisson process 
whose rate depends on the external signal s(t) . The transi- 
tions from 2 to 3 and 3 to 1 are governed by arbitrary WTDs 
W2(t) and W3(t), which do not depend on the external signal. 

The state of our model is described by the three con- 
ditioned probabilities Pi(t\to), i = 1, 2, 3 that the system 
is in state i at time t given that it has been in state 1 at 
time to. These probabilities obey the generalized master 
equations 

A(*|*o) = /dT7(T)PL(r|*o)(tfia ° Ws)(t -r) - 7(*)i\(t|t ) 
A(*|*o) = 7(*)*M*l*o) - [dT 7 {T)Pi(T\t )w 2 (t - t) (1) 



Ps(t\t ) = /dT 7 (r)Pi(T|to) [w 2 (t - r) - (w 2 o 103) (i - r)] 



with initial condition Pi(*o|*o) = 1, jFb(*o|*o) = 0, 
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-P3(*o|io) = and (w 2 ° w 3 )(t) := f Q dTw 2 (T)w 3 (t - r). 
The meaning of eqs.Q will be exemplarily explained by 
means of the first one: The change in probability to be in 
state 1 is the probability influx minus the probability out- 
flux. The later is the instantaneous rate j(t) times the 
occupation probability P\(t\to), as the transition from 
state 1 to state 2 is a Poisson process with rate 7(t). The 
influx into state 1 is related to the outflux of state 1 in the 
past: Consider a system which has left state 1 at some 
instant t in the past, i.e. between to, when the system 
has started in state 1 and the current time t. If this sys- 
tem waits the time t — r in state 2 and 3, which happens 
with probability (w2°Ws)(t— r)), it will reenter state 1 at 
time t, therefore contributing to the influx into state 1 at 
time t. Taken together this leads to the presented influx 
term. The equations for P 2 (i|io) and P3(t\to) follow from 
similar considerations. 

In the case of a constant rate those equations can be 
shown to be equivalent to so called semi Markovian mas- 
ter equations used to describe continuous time random 
walks [1J|. 

For the time dependent rate we use a Kramers type 
rate, which is modulated by the external signal, j(t) = 
ro(D) exp(— cosfH), where D is the noise strength, 
ro(-D) oc exp(— AJ7 e ff/D) is the rate without driving and 
A and f2 are the effective amplitude and frequency of 
the driving. 

In order to show the effect of SR in our model we 
calculate the SNR and the SPA for small a :— Aq/D. 
We assign a high output x(t) = x\ to the system if 
it is in state 2 (firing), and a low output x(t) = .t 
if it is in state l(rest) and 3(refractory) respectively. 
The spectral power of the output x(t) at the input fre- 
quency fl can be derived from the asymptotic mean value 



(x(t)) 



asy 



xo(Pi{t)a S y + P3{t) aS y) + xiP2(t) asy where the 



asymptotic probabilities Pi(t) asy 
are governed by 



lim 



t - 



Pi(t\t ). 



A(*W = -7(*)iM*W+ (2) 

oo 

dry(t - r)P x {t - T) asy (w 2 o w 3 )(t) 



together with the normalization condition Pi(i) aS y 
P2{t) asy + PsWasy = 1 where 

P2(t) asy = / dT 7 (<-T)Pl(t-T) asy Z 2 (r) 

Jo 



(3) 



asy — 

OO (• oc 



(4) 



o Jo 



dr dT'j{t -T- r')Pi(* -T- t')^W 2 {t')z 3 {t) 



Zi(r) := 1 — f Q r dT'wi(r') denotes the probability to wait 
longer than r in state i. Linearizing the rate j(t) with 
respect to a, 7(i) = ro(D)(l — acos(ttt)) and making the 
Ansatz Pi(t) asy — J2T=-oo P k exp(z/cfit) eventually leads 
to a tridiagonal recurrence relation for the Fourier coef- 
ficients pk- Solving this tridiagonal recurrence relation 



following we eventually arrive at 



Pa 

Pi =P-i 
with 



w-'i 



0{a 2 



.a(l-wiffi))(l-w 2 (»)w 3 (»)) 
'2 VLw(l - wi(Q)w) 2 (fi)u)3(Q)) 



(5) 
0(a 3 ), 



u>i(Q) = / dre 1 T Wi(r) and Wi — drrwi(r). (6) 
Jo Jo 

where w = Wi + w 2 + W3 is the mean time for one cy- 
cle of the undriven system and u>i and vbi (fl) are calcu- 
lated according to eq. © with the unperturbed WTD 
for state 1, wi(r) = 7oexp(— 7ot). Higher Fourier coef- 
ficients pk, k > 2 are at least of order 0(a 2 ). Therefore 
one arrives at 



(x(t)) 



asy 



(x) Q + a 



(Xi - X ) 

Qui 



Acos(nt + 0)) + O{a 2 ) 



where (x) Q — 4(ieo(u>i + W3) + 11W2) is the stationary 
solution of the process without external signal and A = 
|r(0)|, $ = arg(r(0)) and 

[ ' 1 - wi{Q)w 2 (n)w3(n) ' [ ' 

The SPA 77 is the ratio between the power of the output 
x{t) at the frequency of the input signal and the power 
of the input signal. It is now given by 

(8) 



V 



(xi 



xo) 



|r(0)| + 0(a 2 



Note that not only w depend on the noise strength D but 
also T(fl) due to the noise dependence of wi(fi). 

The SNR is defined as the ratio between the spectral 
power at the driving frequency and the spectral power 
density in the neighborhood of the driving frequency. 
In lowest order in a this is equal to the ratio between 
the spectral power of the driven process and the spec- 
tral power density of the undriven process, both taken at 
the driving frequency. The spectral power density of the 
process without signal can be calculated using a formula 
from renewal theory |l d | 



Sb(fi) 



4(xi -x ) 2 ReG(ft) 
ID ^ 



with 



G(n) = 



(i-u>a(n))(i-ti>i(n)ti>8(n)) 



1 - wi(fl)w 2 (fl)w 3 (£T) 
Eventually the SNR reads 



SNR = ^o^y + 0(« 4 ) 

1 - U>!(fi) 



(9) 



(10) 



(11) 



4u> 



\G(n)\ 2 



+ 0(a 4 ). 
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For a fixed firing and refractory time, which is a good 
approximation for an excitable stochastic FitzHugh- 
Nagumo system in the low noise regime we have w 2 (r) = 
6(t — Tf) and w 3 (t) = 8(t — T r ) and therefore u) 2 (fi) = 
cxp(— iflTf) and W3(fl) = cxp(— iQ,T r ). This leads to a 
signal to noise ratio which is independent of the driving 
frequency, 



SNR = 



2w 



(12) 



In this case, the SPA is proportional to the spectrum of 
the unperturbed process, 



1 (Xi - X ) 2 

D 2 w 2 



(13) 



2 sin 



2 nTj_ 



2 



+ 7o C 1 - cos(fi(I> + T r ))) + 7o ft sin(fi(T / + T r )) 



with proportionality constant (2D 2 w 2 )~ 1 . Taking into 
account the dependence of 70 and therefore of w on D, 
this function has a maximum at a finite value of D (see 
Fig. [21 which is characteristic for SR. However also as a 
function of the driving frequency the SPA shows maxima 
at certain frequencies, which has also been observed in 
bistable systems with delay Q. In the Markovian case, 



Spectral power amplification 



a 




Figure 2: Spectral power Amplification for fixed firing and 
refractory times, horizontal axes: log 10 D vertical axes: 
Other Parameters: 2> = 1, T r = 2, AU = 0.02, x\ = 
1, xo — 0; The minima, where the SNR vanishes, are located 
at fl = 2nn/Tf, n £ N The maxima are approximately lo- 
cated at n = 2im/(Tf + T r ). 



This SNR depends on the driving frequency Q in a mono- 
tone way. However the SPA 



V = 



(x 1 -x ) 2 (7 3 2 + ^ 2 ) 



D 2 w 2 [{ l0l2 + 7073 + 7372) 2 + (7o + li + iD^ 2 + ^ 4 



can show a maximum with respect to the driving fre- 
quency. A reduction to two states, neglecting state 3 
by choosing 103(1") = S(t) gives a frequency independent 
SNR = ^ and a SPA V = (x^xo) 2 /(D 2 w 2 [( l0 +j 2 ) 2 + 
fi 2 ] which depends monotonously on the frequency. For 
72 = 70 this corresponds to the situation of SR in a 
bistable potential, except that the rate for the transi- 
tion back from state 2 to state 1 is not modulated by 
the external signal, which leads to an additional factor 

1/4 in the SNR compared to the well known expression 

2 

SNR = -^70 for the signal to noise ratio in a two state 
model of a bistable system [l^. 

Another quantity, which shows the influence of noise 
and external driving in an excitable system is the inter- 
spike interval distribution (ISID) W(T). In the case 
of bistable stochastic systems an analogous quantity, 
namely the distribution of times for a transition from 
one stable state to the other and back again, have been 
calculated in ^3]> A n inter-spike interval is de- 
fined as the time between two subsequent firing events, 
which in our model corresponds to two subsequent events 
of entering state 2. We calculate this distribution in the 
limiting case A /D — > 0. The ISID can be written as 



W(T) = / #W(T|0)K0). 



(15) 



where W{T\(j)) is the ISID conditioned by the phase of 
the external signal at the time of entering state 2 and 
p((f>) is the distribution of signal phases at the moment 
the system enters state 2. 

The conditioned ISID W(T\<j>) is given by 



W(T[ 



dt(w 2 o w 3 )(t)wi{T - t\(j> + tn),(l6) 



where wi(t\4>) is the WTD in state 1 given that the phase 
at the moment of entering this state was <j>. This dis- 
tribution reads w\{t\4>) = 70(r)exp(— J Q dtj^,(t)) with 
7^(r) = r o (L>)cxp(-acos(fh- + 0)). 

The phase distribution p(<fi) can be readily calculated 
from the asymptotic solution Pi(t) asy . It is proportional 
to the probability current from 1 to 2 at the correspond- 
ing time t = (f>/Q, which is given by f(t) — 7(£)Pi(t) asy . 
The proportionality constant is fixed by normalization, 



where all transitions are rate processes, i.e. w 2 {t) — 
72exp(-7 2 r) and w${t) — 73exp(— 73T), we get 



2tt 



ll W 



0(a 2 ). (17) 



SNR = 



7l 



n 2 



4w 7^ + 7073 + 73 + ^ 2 



(14) 



The probability current f(t) integrated over one period of 
the external signal gives the average number of entrances 
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into state 2 during one period of the external signal |2fl |. 
As a necessary condition for frequency synchronization 
with the external signal this number has to be one, which 
leads to the time scale matching condition w = . 

Completing the calculation of £>(</>), we insert Pi(t) asy 
using eqs. (01 into eq. I|17fl and arrive at 

pU) = -aAcos((f> + 4>)) +0{a 2 ). (18) 

with A = \g(d)\, 4> — arg(g(fi)) and 



i - wi(n) 



(19) 



The first order in a of p((j>) is sufficient to calculate the 
WTD W(T) according to eq. <dJ| up to order 0{a 2 ). 
This can be seen by writing down the general form of 
p(<t>) up to order 0(a) 2 , p(cj)) = ± + a/i(0) + a 2 / 2 (0). 
Due to the normalization the integral over </> from to 
2-7T of /i and f 2 must vanish. Inspecting Eq. 1151 we 
notice that the term of order 0(a 2 ) which stems from 
0{a 2 ) of p((f>) and 0(a°) of W(7» vanishes as W(T\(f>) 
does not depend on (f> in O(a ). Inserting eq. Ill (it and 
eqs. I|18|l into 11 "il one can at least in principle calculate 
the inter-spike interval distribution for arbitrary WTD 
u>2 and W3. 

As an example we explicitly calculate the ISID for fixed 
waiting times Tf and T r in state 2 and 3 respectively 
which is 



r e 



-r (T-T f -T r )(l+z r ) 



W{T) 

2(l + /3 2 )(l + cos(ftT + 0)) +0{a A ) 



1 



a 2 . 

1 + T^i-2/ 



(20) 



if T > T r +T f and otherwise, tan0 = 2fg/{f 2 -g 2 ), f = 
l+/3 2 (l-cosA0)+/3sinA0and3 = f3 cos A0+/3 2 sinA</>, 
Acj) = Q(T r + Tf). A plot of the resulting ISID is shown 
in Fig. 03 The probability that the interval between 
two subsequent spikes is less than Tf + T r is zero. The 
distribution for larger intervals is an exponential decay 
modulated by the external signal, which leads to small 
peaks at T — 27rn/Q. As the calculation is limited to 
small a, which corresponds to a weak modulation of the 
transition rate from 1 to 2 due to the external signal, 
those peaks are poorly pronounced. 

In conclusion we have derived analytical expression 
for the spectral power amplification and the signal to 
noise ratio for small signal amplitudes in a non Marko- 
vian three state model for stochastic excitable systems. 
The SNR and SPA show a maximum for a certain finite 
value of the noise, which indicates SR. However, for many 
choices of the WTDs in state 2 and 3, the SPA as well as 
the SNR show one or several maxima also as a function 
of the driving frequency, i.e. a "bona fide" resonance. In 
the special case of fixed waiting times in state 2 and 3 the 
SNR does not depend on the frequency, indicating, that 




Figure 3: Inter spike interval distribution W(T) as a function 
of the driving frequency fl. Tf = 1, T T = 2, ro(D) = 1, a = 
0.5 



all signals are equally well separated from the noisy back- 
ground, independent of their frequency. We have further 
derived expressions for the inter-spike interval distribu- 
tion for small signal amplitudes. 
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